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Field  emitters  rely  on  high  gradients  on  microfabricated  structures  to  enable  substantial  levels  of 
emitted  current,  which  can  in  turn  act  to  reduce  the  held  at  the  emission  site  of  a  single  emitter.  An 
account  of  that  effect  is  obtained  in  the  following  steps:  a  model  of  the  emitter  geometry  based  on 
the  point  charge  model  that  allows  for  the  determination  of  the  apex  radius  and  held  enhancement 
factor  for  arbitrarily  sharp  emitter  structures  is  given,  followed  by  an  analytical  formula  for  the 
calculation  of  total  current  from  such  a  structure  and  then  by  a  model  of  the  effect  of  emitted  current 
in  suppressing  the  held  at  the  emission  site.  Predictions  of  the  impact  of  space  charge  on  the  emitted 
current  are  made  and  compared  to  hndings  of  Barbour  et  al.  [Phys.  Rev.  92,  45  (1953)]  for 
sharpened  structures  with  varying  work  function.  A  discussion  is  given  of  the  method  to  combine 
single  tip  three-dimensional  results  with  a  study  of  space  charge  on  held  emission  in  one  dimension, 
[doi:  10. 1063/1.3272688] 


INTRODUCTION 

Field  emitters  are  high  brightness  electron  sources  under 
investigation  for  a  variety  of  applications,  such  as  electron 
beam  lithography  ’  and  transmission  electron  microscopes/ 
spacecraft  propulsion,4'5  millimeter-wave  vacuum  electronic 
amplihers  and  terahertz  devices,6-*  and  advanced  particle  ac¬ 
celerators  and  free  electron  lasers  (FELs).9-13  Their  operation 
can  be  affected  by  the  helds  from  the  charges  they  emit  if 
those  helds  reduce  the  extraction  held  at  the  emission  sites, 
and  such  consequences  are  termed  “space  charge  effects.” 
Such  effects  can  induce  undesirable  growth  in  the  dimen¬ 
sions  of  the  beam  as  it  propagates,  as  commonly  accounted 
for  in  the  beam  envelope  equation;14  a  related  metric  of  beam 
quality  in  that  equation  is  called  “emittance”  and  refers  to  the 
tendency  of  a  beam  to  spread  as  it  propagates  due  to  trans¬ 
verse  velocity  components  in  the  electron  distribution  func¬ 
tion.  An  account  of  both  space  charge  and  emittance  from 
held  emitters  is  required  to  model  how  the  beam  changes  as 
it  propagates.  In  the  present  study,  an  examination  is  under¬ 
taken  of  the  space  charge  forces  that  can  exist  in  the  vicinity 
of  a  single  emitter;  separate  studies  consider  space  charge 
effects  associated  with  an  array  of  emitters15  and  estimations 
of  the  emittance  of  a  single  tip  and  an  array.16 

In  vacuum  electronic  devices,  particle  accelerators,  and 
free  electron  lasers,6-*  electron  sources  are  often  run  space 
charge  limited  (in  the  case  of  thermal  sources)  or  at  such 
high  current  densities  (in  the  case  of  photocathodes  in  rf 
injectors17-19)  that  an  understanding  of  space  charge  and 
emittance  growth  at  the  cathode  is  essential.  As  increasingly 
brighter  beams,  shorter  bunches,  and/or  emission  modulation 
are  sought,  field  emitter  arrays  have  been  under  consider¬ 
ation.  The  simulation  of  intense  sources  is  the  domain  of 
highly  capable  particle-in-cell  (PIC)  simulation  codes20-23 
which  do  have  field  emission  algorithms”4  but  in  the  case  of 
field  emission  from  realistic  sources,  accounting  for  space 
charge  is  hampered  by  the  many  orders  of  magnitude  dispar- 
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ity  between  the  dimensions  of  the  emitting  structures  and  the 
beam  acceleration  and  transport  region  [not  to  mention  the 
inherently  three  dimensional  (3D)  nature  of  an  array  of  field 
emitters  is  in  contrast  to  the  usual  two  dimensional  (2D) 
rotational  symmetry  of  the  electron  beam].  While  in  prin¬ 
ciple,  meshing  a  grid  to  the  dimensions  of  the  emitter  site  is 
possible,  the  current  density,  grid  and  macroparticle  size,  and 
time  step  (if  non-static)  are  inter-related  such  that  computa¬ 
tional  overhead  demands  make  fine  grids  a  burden,  and  so  a 
one  dimensional  (ID)  Fowler-Nordheim  (FN)  equation  uti¬ 
lized  over  the  area  of  the  cathode  with  a  factor  relating  sur¬ 
face  field  to  anode  voltage  is  often  employed.  Even  though 
the  promise  of  high  brightness  electron  sources  such  as  field 
emitters  generates  a  strong  incentive  to  model  them  in  beam 
codes,  the  magnitude  of  the  current  density  at  the  emission 
site  indicates  ignoring  space  charge  effects  on  emission  and 
emittance  for  operational  performance  levels  is  ill  advised. 

A  coherent  theory  of  field  emission  topical  to  vacuum 
electronics,  accelerators,  and  free  electron  lasers  in  particular 
requires,  first,  an  account  of  space  charge  and  emittance  for 
such  sources  under  macroscopic  and  possibly  low  dimen¬ 
sional  environments,  second,  an  account  of  the  microscopic 
and  fully  3D  regions  near  the  emitter  site,  and,  third,  a  char¬ 
acterization  of  the  beam  that  is  to  be  produced  (emittance). 
The  present  work  addresses  the  second  requirement,  as  the 
first  and  third  are  the  subject  of  separate  studies,15’16  al¬ 
though  a  consideration  of  the  methods  to  hand  off  electron 
distributions  from  the  microregime  (the  unit  cell  or  single 
emitter  region)  to  the  macroregime  (the  anode -cathode  or 
beam  transport  region)  is  common  to  all. 

Therefore,  here,  a  methodology  for  the  analysis  of  space 
charge  forces  on  the  emission  process  within  the  unit  cell  is 
described,  and  a  manner  in  which  ID  methods  may  be 
brought  to  bear  on  arrays  of  emitters  operating  together  is 
presented.  Apart  from  providing  a  comprehensive  theoretical 
account,  the  approach  presented  is  intended  to  provide  a 
framework  to  investigate  field  emitters  without  intensive  nu¬ 
merical  efforts  in  a  manner  amenable  to  PIC  codes  for  when 
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space  charge  is  an  issue  and  the  cathode  area  is  but  a  small 
region  of  the  simulation.  The  theory  unfolds  as  follows.  First, 
a  point  charge  model  (PCM)  is  modified  to  treat  a  well- 
studied  class  of  emitters  otherwise  known  as  Spindt-type 
field  emitter  arrays.”  The  revised  PCM  overcomes  the 
disparity  in  length  scales  between  the  emission  site  and  the 
overall  dimensions  of  the  emitter  to  predict  apex  field  and 
field  fall-off  quite  accurately.  Next,  an  analytical  method  to 
calculate  the  total  emitted  current  from  a  single  tip  is  devel¬ 
oped.  From  the  total  emitted  current,  an  estimate  of  the  im¬ 
pact  of  space  charge  on  a  single  emitter  is  developed.  The 
“macroscopic”  field  near  the  cathode  boundary15  functions  as 
the  “far  field”  in  the  vastly  smaller  unit  cells  under  consid¬ 
eration  in  the  present  study.  Finally,  a  comparison  of  the 
single  emitter  space  charge  effects  is  made  with  a  study  by 
Barbour  et  al.2>  in  which  changes  in  I(V)  characteristics 
changed  as  the  work  function  of  the  emitter  apex  was  sys¬ 
tematically  altered  through  the  application  of  barium. 


RELATION  TO  THE  ID  THEORY 


The  potential  away  from  an  array  of  emitters  (ungated  or 
otherwise)  rapidly  becomes  uniform  at  scales  larger  than  the 
tip-to-tip  spacing.15  The  nature  in  which  the  field  near  the 
emitter  site  changes  as  the  distance  from  the  emitter  apex 
becomes  comparable  to  the  tip-to-tip  spacing  is  therefore 
topical.  In  the  PCM  emitters  are  approximated  as  aggregates 
of  point  charges.  Presently,  the  discussion  here  needs  only  to 
consider  a  single  point  charge  to  advance  the  argument,  and 
so  the  question  becomes  how  the  field  varies  for  a  single 
point  charge  as  compared  to  an  array.  The  constant  back¬ 
ground  field  F  of  the  ID  theory  is  intuitively  related  to  how 
the  field  changes  as  the  distance  from  the  array  of  emitters 
increases,  and  showing  that  relation  proceeds  as  follows.  In 
units  of  the  tip-to-tip  spacing  l,  the  width  of  the  array  is  M, 
the  anode  to  cathode  separation  is  N/2,  and  it  is  assumed  that 
N<M  to  approximate  the  ID  configuration.  The  anode  is  at 
Nl  2  to  make  the  potential  there  explicitly  nonvarying  in  the 
transverse  coordinates  and  is  accomplished  by  charges  (es¬ 
sentially  image  charges)  at  N.  The  field  between  oppositely 
charged  planes  of  point  particles  is  then 


z 


v'm2-/ 

k=-\jM2-j2 

( N-z )  j 

2 \M2  +  (N—  z)2  J 


(1) 


where  the  top  line  is  the  sum  over  discrete  point  charges,  the 
bottom  line  is  the  on-axis  large  M  limit  in  which  the  sums 
are  replaced  by  integrals,  NI2<^M  is  the  cathode  to  virtual 
anode  separation,  X  is  the  magnitude  of  the  individual  point 
charges,  F0=kl Atte012  is  the  infinite  plane  (capacitor  plate) 
field,  <p  [with  the  coefficient  of  the  summation  in  Eq.  (1)]  is 
proportional  to  the  potential  associated  with  an  anode- 
cathode  pair  of  point  charges,  or 


<PiM  = 


1 

\/(/  -  x)2  +  O'  -  y)1  +  z2 


yj(i-x)2+  {j  -y)2+  (N-z)2 

and  c p'ij=dz(pjj  (the  z  component  of  the  field  is  the  largest 
component,  and  therefore  considered  by  itself  for  ease).  The 
remaining  terms  are  defined  in  Table  I.  For  very  large  M, 
augmenting  a  transverse  coordinate  (x  or  y)  by  /  does  not 
measurably  affect  the  potential  (for  M=°°,  it  has  no  effect  at 
all),  so  it  shall  be  assumed  that  x  and  y  are  both  smaller  than 
1 12  in  units  of  /.  It  is  expected  that  near  the  center  emitter  for 
z<  1  (that  is,  within  the  “unit  cell”),  the  z  component  of  the 
field  for  the  center  point  charge  should  dominate,  but  that  for 
z>  1,  then  the  field  should  be  significantly  more  insensitive 
to  variations  with  z  and  less  directly  correlated  with  a  single 
charge.  This  intuition  can  be  shown  valid  by  considering  the 
case  M  =  2N=40,  and,  in  particular,  the  z  component  of  the 
field  of  the  point  charge  and  the  array  of  charges  for  on  axis 
(x=y  =  0)  and  at  the  lowest  point  of  the  potential  off  axis 
(x=y=  1/2).  As  shown  in  Fig.  1,  the  center  point  charge  field 
fades  as  1  lz2  until  at  the  unit  cell  boundary  z  ~  1 ,  it  is  over¬ 
shadowed  by  both  the  on-and  off-axis  fields  due  to  the  array, 
which  have  asymptotically  approached  the  analytical  form 
shown  in  Eq.  (1).  Therefore,  the  unit  cell  is  defined  by  the 
region  -l/2<x,y  ^  1/2  and  0<z^  1,  as  schematically  sug¬ 
gested  in  Fig.  2.  The  z=  1  boundary  can  be  assumed  to  be  a 
constant  over  the  region  of  the  virtual  cathode  which  takes 
its  place,  and  the  virtual  cathode  can  then  serve  as  the  bound¬ 
ary  for  the  ID  studies.15  In  such  a  way,  the  considerations  of 
space  charge  near  the  emitter  (unit  cell  and  3D)  can  be  sepa¬ 
rated  from  the  anode-cathode  (ID)  region  that  is  the  normal 
realm  of  the  simulations  modeling  space  charge  limited  cur¬ 
rent. 


THE  3D  MODEL 

A  limitation  of  a  ID  field  emission  model  in  the  analysis 
of  space  charge  is  that  generic  field  emitters  emit  over  square 
nanometer  scale  sites  that  are  spatially  separated  (e.g..  Table 
III  of  Ref.  30)  because  the  sites  have  a  high  degree  of  cur¬ 
vature  to  enable  fields  of  GV/m  from  background  fields  of 
order  0(100  MV/m),  potential  differences  on  the  order  of 
10  kV  for  ungated  geometries  or  100  V  for  gated  geom¬ 
etries.  The  impact  of  space  charge  on  field  emission  from 
multidimensional  geometries  is  difficult,  not  only  because 
finding  the  field  gradient  along  a  curved  surface  over  which 
only  a  small  portion  emits  is  difficult  but  also  because  avail¬ 
able  emission  equations  are  inherently  ID  by  virtue  of  the 
solution  of  Schrodinger’s  equation  for  current  density.  There¬ 
fore,  the  emitted  current  is  obtained  by  considering  modifi¬ 
cations  to  the  parallel  plane  case  (as  done  by  Barbour  et  al. 29 
and  Lau  et  al.31)  or  by  integrating  the  local  ID  current  den¬ 
sity  for  a  given  field  gradient  over  a  multidimensional  sur- 

32 

face  of  the  conductor. 

Finding  solutions  to  the  field  variation  over  the  surface 
of  the  emitter  is  not  trivial.  Exact  conformal  methods  can  be 

33  34 

used  in  2D,  ’  as  can  prolate  spheroidal  techniques  in 

35—38  39—41 

3D.  Numerical  methods  are  available^  with  effort. 
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TABLE  I.  List  of  symbols  and  parameters. 


Symbol 

Definition 

Value 

Unit 

Fundamental  constants  and  parameters 

m 

Electron  mass 

510  999 

eV/c2 

c 

Speed  of  light 

299.792 

nm/fs 

h 

Planck’s  constant 

0.658  212 

ev  fs 

q 

Unit  charge 

i 

q 

Rydberg  constant 

13.605  7 

eV 

Fine  structure  constant 

1/137.036 

So 

Permitivity  of  free  space 

5.526  35  X  10-2 

q2/eV  nm 

Q 

akhc/ 4 

0.359  991 

eV  nm 

F 

Field  at  cathode 

ev/nm 

V 

Potential 

eV 

Copperlike  parameters 

a 

Chemical  potential 

7 

eV 

kF 

(2  mfi)inlh 

13.554  6 

1/nm 

VF 

Fermi  velocity 

1.569  19 

nm/fs 

o 

Work  function 

4.5 

eV 

FN  field  emission  terms  for  Cu-like  parameters 

A0 

Jm  parameter  [Eq.  (4)] 

3.138  94  X  10-5 

A/eV2 

B 

JpN  parameter  [Eq.  (4)] 

65.207  3 

eV/nm 

K 

Jfs  parameter  [Eq.  (4)] 

0.772  81 

v(y) 

Elliptical  integral  function 

1-y2 

l-iln(y) 

t(y ) 

Elliptical  integral  function 

-f(y„)  =  1.061 

y„ 

e-m 

0.606  531 

PCM  parameters 

r 

PCM  scaling  parameter 

n 

PCM  number  of  charges 

a0 

Radius  of  emitter  base 

/xm 

Pexp 

dimensioned  radial  coordinate 

nm 

^exp 

dimensioned  axial  coordinate 

nm 

i 

tip-to-tip  spacing 

fjun 

L 

Anode-cathode  separation 

nm,  cm 

a 

dimensioned  apex  radius 

nm 

as 

gete  radius 

/xm 

8, 

time  separation  between  emitted  e~ 

q/h  iP 

fs 

Aip 

current  from  one  emitter 

/xA 

7(Ffip) 

current  density  at  emitter  apex 

A /cm2 

FIG.  1.  (Color  online)  Behavior  of  the  z  component  of  field  for  a  point 
charge  and  an  array  of  point  charges  on  (x=y=0)  and  off  (x=v=l/2)  the 
symmetry  axis  in  units  of  tip-to-tip  spacing.  By  z=  1,  the  variations  that  exist 
in  the  field  due  to  the  discrete  emission  sites  is  largely  (although  not  com¬ 
pletely)  diminished,  suggesting  z=  1  to  be  the  boundary  between  the  unit  cell 
region  and  the  anode-cathode  region. 


Here,  it  is  desirable  to  analytically  obtain  the  field  enhance¬ 
ment  factor,  the  field  variation  over  the  apex  of  the  emitter, 
and  the  potential  variation  on  the  axis  of  symmetry  in  a 
geometry  that  decouples  the  emitter  apex  shape  from  other 
features  of  the  physical  model  (e.g.,  for  a  hyperbolic  cone 
emitter,  the  cone  angle  is  dependent  on  the  anode-cathode 
separation4-).  Consequently,  an  analytical  model  of  the  emit- 


FIG.  2.  (Color  online)  Schematic  representation  of  the  unit  cell  and  ID 
regimes  suggested  by  Fig.  1  for  gated  and  ungated  geometries. 
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ter  is  developed  that  generalizes  one  used  to  study  field  emis¬ 
sion  and  dark  current  in  accelerators  and  referred  to  there  as 
the  PCM.43 

The  image  charge  FN  equation 

The  FN  equation,  as  widely  used,  is  of  a  form  traceable 
to  Murphy  and  Good44  and  relies  on  elliptical  integral  func¬ 
tions  u(y)  and  t(y).  The  best  approximation  to  these  func¬ 
tions  has  recently  emerged42'46  and  leads  to  a  different  form 
of  the  FN  equation  than  the  most  commonly  used  version  in 
which  u(y)  is  assumed  to  be  linear  in  y2.  Using  the  Forbes 
approximation  to  the  elliptical  integral  function  v(y)  and  a 
representative  constant  for  f(y)  (as  in  Table  I),  the  general 
form  of  the  FN  equation  is  (albeit  in  a  slightly  different  form 
than  that  given  by  Forbes47) 

/FN(U)=AU2-Kexp(-BAF),  (3) 

where  A,  B,  and  k  are  constant  (or  nearly  so)  and  given  by 

8  Q  1 2m 

K=9hV~®’ 


FIG.  3.  (Color  online)  Comparison  of  the  monopole  versus  dipole  tips  for 
r— 1  and  n  =  6. 


V(p,z)  =  Foa0V„ 


Pexp  ^exp 
ao  a0 


n 


V„(p,z)  =  -z  +  (p-  +  z ) 


2\-l/2 


Al 


[ p2+(z-zj)2y 


(5) 


g  /$2e6\*_  /  d>2e6\* 

16Tr/i<I>f(y0)2\  4 Q  )  \  4Q  ) 


B=  —  \/2;w$3, 
3  h 


(4) 


and  where  B  is  equivalent  to  fiFN  in  Eq.  40  of  Ref.  48,  ya 
=  e~1/2  =  0.6065,  t(ya)=  1  +  (l/6e)=  1.061,  and  other  notations 
follow  Ref.  43  or  are  in  Table  I.  The  constant  ya  is  such  that 
the  tangent  line  to  the  Forbes  approximation  to  u(y)  is  linear 
in  y2  with  a  slope  of  -1,  which  provides  a  connection  to  the 
usual  (metal)  result  in  that  if  k=  0,  then  the  usual  represen¬ 
tation  is  recovered.  For  copperlike  parameters  it  follows  that 
k= 0.772  81.  While  Eq.  (3)  is  not  in  a  form  resembling  the 
FN  equation  as  originally  derived,  it  is  a  consequence  of  the 
forms  used  for  v  (y)  and  f(y). 

As  a  final  observation,  as  discussed  (and  advocated)  by 
Forbes,47  the  form  of  Eq.  (3)  removes  the  privileged  role  of 
so-called  FN  coordinates  in  which  In (JIF2)  is  plotted  against 
1  IF,  and — by  extension — the  same  with  ln(/tip/V2)  with 
1  /  V.  In  point  of  fact,  ln(7FN)  will  appear  linear  as  a  function 
of  1  IF  for  a  variety  of  k  values.  The  pedigree  of  Millikan- 
Lauritsen  plots,  in  which  ln(7)  versus  II F  is  shown,  is  older 
than  FN,  easier  to  use,  more  robust  in  dealing  with  correc¬ 
tions  and  sources  of  voltage  dependence,  and  therefore  shall 
be  how  I{V )  data  shall  be  represented  here. 


where  X0  is  implicitly  set  equal  to  1  and,  in  the  second  line, 
p  and  2  are  dimensionless  cylindrical  coordinates  and  corre¬ 
spond  to  the  ratio  of  the  physical  radial  and  axial  coordinates 
with  the  unit  a0  such  that  V„(p,z )  itself  is  dimensionless.  The 
Zj  are  defined  according  to 


n- 1 


I  -  /  " 


1  —  r 


Sn(r) 


(6) 


and  denote  the  location  of  the  nth  charge  of  relative  magni¬ 
tude  X„.  The  physical  height  of  an  n -charge  emitter  is 
Sn+l(r)aQ.  The  factor  r  is  a  scaling  parameter  such  that  a 
charge  Xj  is  placed  F  above  the  one  below  it  (i.e.,  z,J+ \  =  F 
+ z,j) .  The  values  of  X;  are  dictated  by  the  boundary  condition 
V/„(0,z„+])  =  0,  and  the  surface  is  defined  by  the  zero  equipo- 
tential  line.  Because  charges  of  only  one  sign  are  considered 
in  Eq.  (5),  V„  in  Eq.  (5)  shall  be  termed  the  “monopole” 
model  V™"m'(p,z)  to  distinguish  it  from  another  version  be¬ 
low.  The  monopole  emitters  have  broad  bases  suitable  to 
model  cone  formation  from  heating  and  melting.  In  contrast, 
microfabricated  field  emitters  tend  to  be  of  a  sharper  conical 
shape  and  carbon  nanotubes  more  cylindrical.49  It  is  there¬ 
fore  useful  to  introduce  a  model  that  more  closely  approxi¬ 
mates  narrow  geometries  and  called  the  “dipole”  model 
V^p\p,z)  because  it  considers  pairs  of  opposite  charges 
placed  symmetrically  about  the  z  plane.  The  dipole  Vn  is 
defined  by 


The  numerical  model  for  a  single  emitter 

The  PCM  gives  field  enhancement  and  tip  radii  param¬ 
eters  in  a  computationally  expedient  manner.  A  line  of 
charges  is  placed  along  the  z  axis  above  the  z=0  plane  in 
such  a  manner  as  to  have  the  zero  equipotential  line  approxi¬ 
mate  the  shape  of  an  emitter.  The  potential  everywhere  is 
related  to  the  background  field  Fa  (in  practice,  the  field  at  the 
virtual  anode  in  the  unit  cell)  to  the  base  dimensions  of  the 
emitter  characterized  by  a  length  scale  aQ  via 


vT (pA  =  -  7  +  2  M[p2  Az  -  z/]-1/2 

i=  i 

~lp2  +  (z  +  Zj)2T112},  (7) 

where  X0  is  implicitly  set  equal  to  0.  It  is  subject  to  the  same 
boundary  condition  \4jlpl(0,z„+l)  =  ().  The  presence  of  mirror 
image  charges  across  the  ~=0  plane  results  in  a  much 
sharper,  more  ellipsoidal  emitter  shape,  as  shown  in  Fig.  3 
for  the  particular  case  of  r=  1  and  n  =  6.  In  both  the  monopole 
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and  dipole  models,  the  field  enhancement  factor  f3n  and  apex 
radius  a„  are  given  by 

/3„(r)  =  -  dzVn(0,z)\z=Zn+1, 


a„(r)  = 


<?pVn(p,z) 


P= 0,z=z„+i 


(8) 


It  is  observed  that  Eq.  (8)  refers  to  apex  values  of  field  en¬ 
hancement  and  radius,  but  Eq.  (7)  does  not  preclude  the 
evaluation  of  such  quantities  off-axis  (rather,  the  on-axis  val¬ 
ues  are  the  most  appropriate  metric  to  discuss).  The  field  and, 
hence,  local  field  enhancement  are  readily  evaluated  from 
Eq.  (7)  off  axis  and  down  the  sides  of  the  emitter.  Impor¬ 
tantly,  the  PCM  formulation  allows  for  the  treatment  of  small 
asymmetric  protrusions  existing  on  a  conical  emitter  through 
the  addition  of  an  off-axis  point  charge  in  a  straightforward 
manner  or  alternately  a  whisker  that  is  bent  and  therefore  no 
longer  rotationally  symmetric.  These  generalizations  are  use¬ 
ful  extensions  of  the  theory,  but  are  not  presently  considered. 

The  methods  to  evaluate  /3„,  and  an  in  terms  of  r  and 
n  are  now  given.  As  may  be  anticipated,  \„(r)  for  small  r  is 
an  exponentially  shrinking  quantity  as  n  increases,  and  it  is 
computationally  expedient  to  introduce  an  auxiliary  function 
Pn{r)  defined  by 


Kir)  =  P„(r)r" 


(9) 


to  separate  out  the  exponentially  small  behavior.  The  bound¬ 
ary  conditions  can  now  be  concisely  expressed  by  Eq.  (5) 
and  (7)  evaluated  at  p=0  and  z=zn+\  via 


C 

0 

E 

0 

O 


T3 
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FIG.  4.  (Color  online)  Dimensionless  field  enhancement  factor  f3n  and  apex 
radius  factor  an  as  a  function  of  n  using  the  dipole  model  of  Eqs.  (13)  and 
(15).  The  numbers  following  (5  and  an  in  the  legend  are  the  value  of  r  for 
that  line. 


The  matrix  method  can  likewise  be  used  to  calculate  the 
field  enhancement  factors  (3n(r),  which  like  Pn(r)  can  be  cast 
as  the  solution  to  matrix  equations.  Define  B  by 


r  d  monol  _ 

J  nj- 


1 

W(S„+ w)2  ’ 


[Bdipl]nj-  = 


4SjS  „+i 


^n+i-/5„+i+5J.)]2’ 


(12) 


where  the  argument  in  Sj  is  suppressed.  Then 

^mono  _  gmono  .  p  _|_  Y 


yjdipl  _  gdipl  .  p  +  J 


03) 


^M^no(r)Pj(r)=Sn+l(r)- 


7=1 


s„+iM’ 


^M*f(r)Pj{r)=Sn+1{r). 

7=1 

In  matrix  parlance, 

[Mm0n0]nJ.  : 


0„ 


3'/i+ 1  -jir) 


I  I  .=  -  2Sj(r)  _  q 

L  in’J  Sn+i_j(r)[Sj(r)  +  5„+1(r)] 


00) 


(ID 


where  I  is  the  identity  matrix  and  [Y]„=  1  +[5„(r)]  2.  As  with 
M,  B  is  a  lower  triangular  matrix  and  so  all  the  field  en¬ 
hancement  values  for  any  n  =£  N  can  be  obtained  by  inverting 
the  matrix  equation  for  n=N.  The  equations  for  apex  radius 
can  be  put  in  the  same  form,  but  it  is  convenient  to  simply 
give  them  as,  for  the  monopole  case  (the  r  dependence  of  the 
S  terms  is  suppressed  for  visual  clarity). 


a 


mono 

n 


=  s 


n+ 1 


2 


(14) 


and  for  the  dipole  case 


where  0n  ;  =  1  if  y  =£«  and  0  otherwise,  and  use  has  been  made 
of  Sn+1-Sj  =  PSn+ 1  _j.  Equation  (10)  is  therefore  concisely  ex¬ 
pressed  as  matrix  equation  M  P=S,  where  M  is  a  lower 
triangular  matrix,  and  therefore  the  inversion  P=M  ■'•S  is 
easily  accomplished  using  standard  numerical  techniques. 
Because  M  is  lower  triangular,  the  solution  for  n=N,  where 
N  is  larger  than  any  n  likely  encountered  in  practice  (e.g.,  at 
A=24  and  r=0.4,  the  base  of  a  nanometer-radius  emitter 
would  be  approximately  2  m),  need  be  performed  but  once 
for  a  given  r  and  the  coefficients  can  be  stored  and  therefore 
available  for  any  n<N  in  subsequent  calculations.  Such 
would  be  useful,  for  example,  to  model  the  physical  growth 
of  an  emitter  as  n  increases. 


P[ _ SjSn+i _ 

dipi  7=1  ^  [3„+i-;(S;  +  3„+i)] 

an  ~  =  n  /  9  o 

■y*  P j  Sj(S,  +  3S„+1) 

~i  i-22[3„+i-j(3J  +  5)I+1)]2 


(15) 


The  behavior  of  the  field  enhancement  factor  and  the  apex 
radius  for  various  r  and  n  is  shown  in  Fig.  4.  To  relate  to 
physical  emitters,  the  overall  shape  of  the  emitter  is  gov¬ 
erned  by  type  (cusplike  corresponds  to  “mono,”  conical,  or 
ellipsoidal  to  “dipl”),  height  is  governed  by  Sn(r)a^p,  and 
apex  radius  is  governed  by  an(r)a^p,  where  g!qXP  is  a  dimen¬ 
sioned  unit.  In  contrast  to  all  other  analytical  emitter  models. 
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FIG.  5.  (Color  online)  A  dipole  tip  for  r=0.75  and  increasing  values  of  n, 
demonstrating  how  the  shape  of  the  tip  stays  approximately  the  same  even 
as  the  apex  becomes  sharper  for  increasing  n. 


the  model  retains  the  same  basic  emitter  shape  as  the  tip  is 
progressively  sharpened.  This  is  seen  in  Fig.  5  as  n  increases 
for  a  given  r. 


The  tip  current  model 

Near  the  apex  of  the  emitter,  the  surface  is,  to  a  good 
approximation,  parabolic  and  therefore  defined  by  the  equa¬ 
tion 

2 

<Turl<P)  ~  7-n+ 1  “  /.  •  (16) 

2  a„ 

The  current  density  drops  rapidly  as  p  increases  because  to  a 
good  approximation  ,F~,Ftipcos  P=Ttip[l-(p/a„)2]1/2, 
where  the  cylindrical  coordinate  p  is  related  to  the  spherical 
coordinate  0  by  p=an  sin  6.  For  example,  when  p  is  75%  of 
a,„  then  the  current  density  is  1%  of  its  apex,  or  “tip,”  value. 

TO 

The  total  current  is  obtained  from 

f(Ftip)  =  [  JdCl=  f  2 rip \  1  +  (dpZsmf)2J[F(p)]dp 
J  n  Jo 

~  27r(aoXpa„)2y(Ttip)  f  jcV  1  +  x2J[F(anx)\dx. 

Jo 

(17) 

Using  Eq.  (16)  in  Eq.  (3)  gives 
1  /  7 r\  1/2 

7(^tiP)  =  {Erf[/7(F)i(T)] 

-Erf[.v(F)]}exp(.y(F)2)7(Ftip), 


p(x)  = 


5 B  -  4(2  -  k)x 
2 B  -  2(1  -  k)x ’ 


six)  = 


B  -( 1  -  k)x 
\/2x[3 B  -  2(3  -  k)x\ 


(18) 


where  Erf(x)  is  the  error  function  and  the  tip  radius  is  a 
=  a„aQXp  (and  is  therefore  dimensioned).  For  example,  an 
apex  field  of  8  GV/m,  a  tip  radius  of  a^p=  10  nm,  and  cop¬ 
perlike  parameters,  for  which  p  and  s  are  equal  to  2.0945  and 
3.6134,  respectively,  Eq.  (17)  predicts  that  /tip=49.7  pA, 
whereas  Eq.  (18)  gives  52.7  pA.  The  notional  emission  area 


FIG.  6.  (Color  online)  The  notional  emission  area  g  factor  [ratio  of  tip 
current  7tip  to  the  product  of  the  current  density  on  axis  with  the  area  of  a 
hemisphere  of  radius  equal  to  the  apex  radius  7rcrJ(Ftip)]  calculated  exactly 
(open  circles),  approximately  using  Eq.  (19)  (dashed  line),  and  by  the  qua¬ 
dratic  fit  of  Eq.  (20)  (solid  line). 


Utip/Aip)  is  0.176  44rra2  corresponding  to  55.431  nm2.  The 
field  dependence  of  the  notional  emission  area  is  shown  in 
Fig.  6.  Using  the  r=0.4  line  of  Fig.  4,  this  would  be  achieved 
by  an  emitter  characterized  by  agxp=  1.5345  pm  and  n  =  6  in 
a  background  field  of  222.61  MV/m.  For  comparison,  such  a 
field  is:  a  fraction  of  the  breakdown  fields  associated  with 
various  metals;  less  than  the  400  MV/m  surface  field  that 
characteristic  of  the  compact  linear  collider  and  compa¬ 
rable  to  the  ratio  of  the  gate  potential  (75  V)  to  the  gate 
radius  (0.375  pm)  for  a  Spindt-type  field  emitter  operating  at 
50  mA  for  a  10  tip  array."  Continuing,  the  dimensionless 
axis  potential  V„(0,z)  and  its  gradient  for  z^z„+i  are  shown 
in  Fig.  7  (albeit  for  n= 4  rather  than  6  to  bring  out  the  near¬ 
apex  behavior)  for  a  dipole  geometry:  by  z~z„+ 1~  1  (a  dis¬ 
tance  characteristic  of  the  emitter  base  radius)  the  field  gra¬ 
dient  is  only  28%  larger  than  the  background  field  whereas 
the  field  enhancement  factor  is  15.855. 

The  coefficient  of  J  in  Eq.  (18),  or  rather  g(F) 
=  /tip/(7Tfl2/),  is  to  a  good  approximation  a  quadratic  func¬ 
tion  of  F,  albeit  that  its  evaluation  is  somewhat  cumbersome. 
It  might  be  thought  that  a  useful  approximation  can  be  ob¬ 
tained  by  using  asymptotic  expansions  of  the  error  function 
Erf(x)  ~  1  -  7r-l/2e-r  (x-1  —  2x-3)  to  give 

o 
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FIG.  7.  (Color  online)  The  dimensionless  factor  Vn(p,z)  of  Eq.  (7)  (open 
circles)  evaluated  on  axis  (p= 0)  and  its  gradient  (solid  circles)  with  respect 
to  z-  The  gradient  is  within  28%  of  its  asymptotic  value  by  z=l,  indicating 
that  field  enhancement  is  curtailed  at  a  length  scale  on  the  order  of  the 
emitter  size. 
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g(F)~ 


F(2-s(F )~2) 
B  +  (l-  k)F’ 


(19) 


but  this  approximation  is  poor  when  F  approaches 
10  eV/nm,  as  shown  in  Fig.  6.  For  computational  purposes, 
it  is  better  to  exploit  the  quadratic  behavior  explicitly  and 
take  advantage  of  g(0)  =  0  to  obtain 

g(F )  »  F(ciF  +  c0), 


c,=2g(Fmax)-4g(^), 

c0  =  -^mJ  +  4g(^j,  (20) 

where  F'max=  10  eV/nm.  For  copperlike  parameters,  ct 
=  -0.064  203  and  c2  =  0.286  47  and  gives  the  line  labeled 
“Quadratic”  in  Fig.  6,  entailing  that  g(F )  therefore  need  only 
be  computed  twice  for  a  good  parametric  representation. 


The  space  charge  model 

With  an  estimate  of  7tip,  an  evaluation  of  when  space 
charge  forces  affecting  emission  for  a  single  tip  can  be  made. 
Assume  that  electrons  are  emitted  along  the  z  axis  every  St 
=  q/Itip  with  (in  contrast  to  a  previous  treatment  ’6)  an  initial 
velocity  v  =  h(k)lm~hkFl m  =  vF,  where  h2k2Fl 2m =/jl.  16  As¬ 
sume  that  the  field  that  accelerates  the  emitted  electron  is 
constant.  In  actuality,  the  field  drops  away  from  the  emitter 
apex  as  a  consequence  of  geometry  and  so  electrons  will  not 
be  swept  away  as  quickly  as  assumed.  The  effect  of  electrons 
further  out  (where  the  approximation  degenerates)  is  much 
less  than  those  closer  to  the  emission  site  (where  the  approxi¬ 
mation  is  good)  because  their  influence  fades  as  the  square  of 
their  distance.  Therefore,  the  cumulative  space  charge  effect 
will  be  larger  than  that  predicted  by  the  constant  field  ap¬ 
proximation  which  tends  to  underestimate  it  (i.e.,  the  con¬ 
stant  field  result  is  a  lower  bound  estimate).  The  distance  of 
the  j th  electron  from  the  apex  is  then 

Zj(St)  =  —(jSt)  +  ^-(jSt)1.  (21) 

m  2m 


The  sum  of  the  forces  on  the  apex  of  the  emitter  from  the 
emitted  electrons  is 

2  00 

SFe~Y—'2[zj(St)T2.  (22) 

4-7780  M 

Let  the  image  charge  be  that  of  a  sphere  of  radius  a,  for 
which 


O  00 

q  -y Zj(St)  +  a 
A-tteqCI j=i  [Z]m2  ' 


(23) 


Introducing  the  dimensionless  term  7j(F)  =  2hkFl  FSt  then  the 
total  reduction  in  field  at  the  apex  of  the  emitter  due  to  the 
line  of  emitted  charge  is  the  sum  of  Eqs.  (22)  and  (23),  or 


FIG.  8.  (Color  online)  Behavior  of  R  and  R'  [Eq.  (25)  with  the  approxima¬ 
tions  of  Eqs.  (26)  and  (27)  for  n= 2]. 


2 a  f.hcm  1/  aF  \  ,  , ,  I 

(24) 

The  functions  R  and  R'  (where  the  prime  does  not  indicate 
derivative)  are  defined  by 


Fix)  =  2 
j=  i 


1 

jij  +  x)’ 


R'(x)  =  ^~2 - 


;=i  j~(j  +  x) 


2‘ 


(25) 


Special  cases  are  /?(0)  =  ^(2)  =  7t2/6,  and  /?(1)=1,  and  R'( 0) 
=  £(4)  =  7t4/90  and  R'( l)  =  2£(2)-3,  where  £(n)  is  the  Rie- 
mann  zeta  function.  The  summations  can  be  approximated 
by  equivalent  integrals,  but  the  error  is  largest  for  small  j. 
Therefore,  it  is  most  efficient  computationally  to  sum  the 
first  few  terms  and  approximate  the  remainder  by  integration. 
We  use 


R{x) 


1 


2  n(n+x)  x 


n- 1 


-in  i  +  - USttt1 


7=1  j(j  +  X) 


(26) 


For  n= 2,  Eq.  (26)  is  accurate  to  within  1.2%.  Performing  the 
same  analysis  with  R'{rj)  results  in  the  approximation 


R'{x) 


2  n  (n  +  x)  (2  n  +  x)  +  x2 
2  n2x2(n  +  x)2 


77-1 


+  2 

7=1 


1 

j2ij  +  x)2’ 


(27) 


which  for  n= 2  is  accurate  to  within  0.6%.  Larger  n  rapidly 
results  in  even  greater  accuracy.  The  behaviors  of  R  and  R' 
are  shown  in  Fig.  8  for  values  of  argument  encountered  for 
copperlike  parameters  (observe  that  n  here  is  not  the  n  of  the 
PCM).  The  behavior  of  SF(F)  is  shown  in  Fig.  9  for  a 
=  10  nm  and  $  =  4.5  eV. 

The  field  at  the  apex  of  the  emitter  is  therefore  smaller 
by  SF  than  the  field  that  would  exist  in  the  absence  of  space 
charge.  The  leading  order  estimate  of  single-tip  space  charge 
is  then  that  the  relationship  between  the  field  F  at  the  apex  to 


Author  complimentary  copy.  Redistribution  subject  to  AIP  license  or  copyright,  see  http://jap.aip.org/jap/copyright.jsp 


014905-8  Kevin  L.  Jensen 


J.  Appl.  Phys.  107,  014905  (2009) 


F  [eV/nm] 


FIG.  9.  Behavior  of  SF  for  a  dipole  geometry  with  a=  10  urn.  work 
function=4.5  eV,  and  calculated  using  Eqs.  (24)  with  St=q/Itip. 

the  anode  voltage  must  be  2/3„VI  Nl=F+  SF.  From  Fig.  4, 
exponential  approximations  for  the  dipole  geometry  for  r 
=  0.4  are  given  by 

/3„(r  =  0.4)  ~  3.125<?0'4095", 

an(r  =  0.4)  «  1.390e-0'9067".  (28) 

The  results  of  the  space  charge  effects  evaluation  are  shown 
in  Fig.  10  for  several  values  of  n  for  a  dipole  configuration 
under  the  assumption  that  the  work  function  is  4.5  eV  and 
the  base  radius  a0  is  1  /tm .  The  Fa  used  as  a  background 
field  has  conceptual  difficulties  when  being  related  to  a  gated 
field  emitter,  as  it  will  bear  a  nontrivial  relationship  to  the 
gate  voltage  and  gate  radius,  and  requires  modeling  specifics 
beyond  the  present  treatment  (but  intimated  elsewhere36). 
The  behavior  shown  in  Fig.  10  for  small  large  macroscopic 
field  [given  by  (F+  SF)  /  /3„]  shows  the  impact  of  tip  sharp¬ 
ening  on  the  onset  of  space  charge  effects  and  bears  a  quali¬ 
tative  similarity  to  data  shown  in  Ref.  27  and  36  (albeit  the 
latter  are  shown  on  traditional  FN  plots):  however,  when 
emission  is  from  an  array,  it  is  far  more  likely  that  space 
charge  between  the  anode-gate  planes  causes  the  character¬ 
istic  turnover  than  the  single-tip  effects  of  Fig.  10.  To  model 
further  space  charge  effects  on  arrays,  the  3D  tip  model  must 
be  embedded  in  the  ID  gap  model,  and  such  an  analysis  shall 
be  deferred  to  a  separate  study. 


FIG.  10.  (Color  online)  Tip  current  on  a  Millikan-Lauritsen  plot  evaluated 
using  Eq.  (18)  (where  increasing  n  refers  to  sharper  emitters  as  per  the 
PCM)  and  for  a  work  function  of  4.5  eV  showing  the  effects  of  emitted 
charge  on  the  FN  behavior  of  a  single  (micro)  emitter  tip  as  a  function  of  n 
for  field  enhancement  and  apex  radius  factors  as  given  by  Eq.  (28).  The 
units  are  amperes  for  current  and  eV/nm  for  field. 


FIG.  11.  (Color  online)  Comparison  of  the  onset  of  space  charge  effects  on 
the  experimental  data  of  Barbour  el  al.  for  a  tungsten  tip  progressively 
covered  with  barium  (l=no  barium;  4=maximum  coverage  of  barium)  as 
shown  on  a  Millikan-Lauritsen  plot.  Barium  causes  the  work  function  to 
decrease,  thereby  changing  the  FN  characteristics. 

Comparison  to  experiment 

Therefore,  for  comparison,  consider  a  sharpened  needle 
rather  than  a  microtip,  and,  in  particular,  for  field  emission 
data  from  a  tungsten  emitter  examined  by  Barbour  et  al. 
(Fig.  3  in  Ref.  29)  in  which  current  intensified  as  barium  was 
added  to  the  needle  to  lower  the  work  function.  The  experi¬ 
mental  data  allow  for  a  comparison  to  the  theoretical  model 
after  the  specification  of  work  function,  tip  radius,  and  field 
enhancement  factor.  First,  the  emission  area  of  the  needle  is 
far  larger  than  for  the  Spindt-type  emitters:  as  a  result,  the 
reported  apex  radius  of  the  emitter  can  be  measured  instead 
of  inferred.  Second,  as  the  uncoated  tip  was  tungsten,  the 
field  enhancement  factor  for  the  needle  can  be  obtained  from 
that  tip  and  used  for  the  coated  tip  measurements.  Third,  with 
the  field  enhancement  factor  and  the  apex  radius,  the  work 
function  for  the  remaining  three  lines  corresponding  to  in¬ 
creased  barium  coverage  of  the  tip  can  be  determined.  It 
needs  to  be  emphasized  that  the  single  unknown  quantity  for 
each  line  (field  enhancement  factor  for  uncoated  needle, 
work  function  for  partially  coated  needle)  can  be  found  from 
a  single  I(V)  data  point  using  Eq.  (18)  rather  than  from  least 
squares  fitting  to  find  slopes  and  intercepts  on  a  FN  plot. 

The  apex  radius  of  the  tungsten  emitter  is  suggested  by 
the  measurements  of  Ref.  29  to  be  364  nm.  Using  an  I(V) 
data  point  extracted  from  line  1  (uncoated  tungsten)  from 
Fig.  3  of  Ref.  29  (V=  7980.6  V,  7=0.646  /jlA),  the  field  en¬ 
hancement  is  determined  to  be  [3=45 1 6q/cm.  Using  single 
data  points  from  the  other  lines  2-4  at  comparable  current 
levels,  and  taking  vF=cl  152.4  (the  chemical  potential  of 
tungsten  is  greater  than  that  of  copper),  the  work  functions  of 
the  partially  coated  tips  are  determined  to  be  3.30,  2.89,  and 
2.52  eV,  respectively.  I(V)  data  are  then  generated  using 
Eqs.  (18)  and  (24).  The  anode  voltage  is  determined  by  V 
=  (F+  SF)  /  [3.  and  the  theoretical  estimates  of  space  charge 
affected  emission  from  the  tip  are  shown,  along  with  the  data 
of  Barbour  et  al.  in  Fig.  11,  where,  following  the  recommen¬ 
dation  of  Forbes,47  a  Millikan-Lauritsen,  rather  than  FN,  co¬ 
ordinate  plot  is  used:  the  numbering  of  the  lines  corresponds 
to  the  same  numbering  as  used  by  Barbour  et  al.  As  noted 
previously,  the  assumption  of  a  constant  (versus  a  declining) 
field  away  from  the  emitter  apex  entails  that  the  magnitude 
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of  SF  will  be  underestimated,  but  it  is  seen  that  even  for 
unoptimized  parameters,  the  prediction  of  the  onset  of  space 
charge  effects  is  quite  clearly  predicted. 

CONCLUSION 

Field  emission  is  particularly  subject  to  space  charge  ef¬ 
fects  because  of  the  strong  variation  in  the  emitted  current 
with  the  field  that  exists.  Because  of  the  high  current  densi¬ 
ties  that  are  possible  with  field  emission,  space  charge  effects 
can  manifest  themselves  in  two  ways:  first,  it  can  suppress 
the  macroscopic  field  between  cathode  and  anode  as  would 
exist  for  an  array  of  emitters  (and  as  described  by  the  usual 
Child’s  law  formulation),  and  second,  for  a  single  emitter, 
charge  near  the  emitter  site  can  reduce  the  apex  field.  In  the 
present  work,  a  theory  for  the  latter  effect  was  developed  and 
used  to  analyze  single  emitter  data  in  the  literature.  It  was 
determined  that  space  charge  can  have  an  impact  even  for  a 
single  field  emitter,  and  a  method  to  estimate  that  impact  was 
given.  Finally,  a  discussion  of  the  method  to  incorporate  such 
a  unit  cell  model  into  a  larger  macroscopic  model,  in  a  man¬ 
ner  suggested  to  be  useful  to  the  needs  of  PIC  codes,  was 
given. 
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